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A simple technique is proposed for numerically determining equilibrium ion distribution functions 
belonging to free energies of the Poisson-Boltzmann type. The central idea is to perform a conven- 
tional Monte-Carlo simulation using the free energy as the "Hamiltonian" entering the Metropolis 
criterion and the spatially discretized density as degrees of freedom. This approach is complementary 
to the possibility of numerically solving the differential equations corresponding to the variational 
problem, but it is much easier to implement and to generalize. Its utility is demonstrated in two 
examples: valence mixtures and hard core interactions of ions surrounding a charged rod. 
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I. INTRODUCTION 

Understanding the interaction between charged 
macroions in solution requires knowledge of how the ions 
from the surrounding electrolyte screen the niacroionic 
charge. A common starting point for the description of 
such systems is a mean field theory, which expresses the 
free energy as a functional of the ionic density. Account- 
ing only for the electrostatic energy and the entropy of 
this density field yields (after functional minimization) 
the nonlinear Poisson-Boltzmann (PB) equation, which 
can be solved exactly in certain special casesETa. Since 
this approach neglects statistical correlations as well as 
the finite size of the small ions, various methods have 
been suggested which approximately account for those 
phenomena by adding correction terms to the free energy 
functionalElEl. However, the analytical treatment of those 
extended Poisson-Boltzmann theories is even harder than 
that of the original one; consequently, the desired ionic 
density profile is computed by numerically solving the 
complicated differential- or even integro-differential equa- 
tions. 

In this paper a method is proposed for determining the 
equilibrium ion distribution, which completely avoids dif- 
ferential equations. The key idea is to perform a straight- 
forward Monte-Carlo (MC) simulation of the free energy 
functional. Although this is a numerical approach as well, 
it has the big advantage that its computational effort 
is mostly unaffected by the mathematical complexity of 
the functional under investigation. Only recently similar 
techniques have successfully been applied to the |Onsager 
free energy functional from liquid-crystal theorja. 



II. DESCRIPTION OF THE MONTE-CARLO 
METHOD 

This section briefly outlines the physical and mathemati- 
cal background and proceeds with describing the Monte- 




FIG. 1. Geometry of the cell model. An infinitely long 
cylinder of radius ro and line charge density A > is coaxi- 
ally enclosed in a cylindrical cell of radius R. Global charge 
neutrality of the system is ensured by adding an appropriate 
amount of oppositely charged counterions. The outer radius 
R is used to match the polyelectrolyte density of a solution 
containing many such macromolecules. 



Carlo scheme. For the sake of concreteness the deriva- 
tions are presented within the framework of the cylin- 
drical cell model, which is a commonly used approach 
for reducing the complicated many body problem of a 
solution of stiff polyelectrolytes (likf, e.g., DNA) to an 
effective one polyelectrolyte theorya. Upon (i) assum- 
ing that the surface of zero electric field surrounding one 
such molecule is on average a cylinder and (ii) neglecting 
edge effects by taking the length of the molecule to be 
infinite, one arrives at the geometry depicted in Fig. |l|. A 
convenient measure for electrostatic interaction strength 
is the Bjerrum length Ib = iSeQ/Aire {k^T = (3~^ being 
the thermal energy and e the dielectric constant), which 
permits the Coulomb energy for two unit charges eo a 



distance d apart to be written as k-oT (.-q/cI. Within a 
mean field density functional approach the counterions 
surrounding the rod are replaced by a cylindrically sym- 
metric ion density n{r) [r is the radial coordinate) which 
gives rise to a likewise symmetric electrostatic potential 
'0(r). The simplest ansatz for the free energy restricts to 
the electrostatic part of the energy and the entropy: 



F[n{r)] ^L dr2Trr 
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Here, n is the average counterion density and L is the 
length of a cylinder subsegment, with respect to which 
all extensive quantities are to be measured in the follow- 
ing. Functional minimization of F under the constraint 
of global charge neutrality leads to the nonlinear Poisson- 
Boltzmann equation 

AV-Cr) = V" (r) + -ij' (r) = n(i?)e''"°'^('") , (2) 

r 

which has to be solved subject to the boundary condi- 
tions 



, i^'iR) = and MR) = 0. 

27rero (3) 
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The analytical solution of (H,H) and in particular a dis- 
cussion of the phenomenon of Manning condensatioix^, 
which occurs for a charge parameter ^ = XI-q/cq larger 
than 1 can e.g. be found in Refs.DU. The present paper 
suggests a method for (numerically) finding the distri- 
bution n{r) which minimizes functionals like (ll) without 
using their corresponding differential equations. 

Imagine the cell being subdivided by the M -f 1 (not 
necessarily equidistant) splitting points sq,si,...,sm 
(with So = ^0 and sm = R) into M concentric cylin- 
drical shells of volume Vi — 7r(s| — s^_i)L. If further- 
more the Ni ions within shell i are replaced by a density 
rii = Ni/Vi, the set {Ni} already completely specifies the 
(macro-)state of the system. The total electrostatic en- 
ergy E of this configuration can be computed exactly by 
piecewise integration of the Poisson equation in cylindri- 
cal geometry, yielding 

M . 
E{m}) = -—Y,\{Q^-l~eon^7Tst^LYln^^ 

I— 1 ^ 

+ eoN,i^Q,^i+eon,T:^ — _^J-LJ I (4) 

where Qi^i is the total charge up to but not including 
shell i. The idea now is to perform a Markov process 
by randomly moving particles between the shells and aer 
cepting the moves with the usual Metropolis criteriontil. 
In the canonical ensemble the probability of a state {Ni} 
is proportional to the product of the Boltzmann factor 
Pe = exp{—f]E} and the probability Ps of distributing 



the (indistinguishable) ions between the shells in a way 
compatible with the given state, i.e.. 



Mm) = ^n^ 

i—1 * 



(5) 



where N = J^i Ni and V = J^i ^- Generate a new state 
{iVj'} from the old one by moving one ion from shell k to 
shell I. Since 



Ps{{Ni}) _ Nk Vi 



,old 



Ps{{N,}) VkNi + 1 nf 



(6) 



detailed balance between the total probabilities finally 
yields the following acceptance probability of such a 



Prob( k^l)= min |l, ^g"''!^"™-^"") 1 



(7) 



There is an alternative way of looking at this: Since the 
entropy of a state is given by 5 = fee In Ps , the expres- 
sion entering in (M) can actually be written as exp{—(3F}, 
i.e., the "combinatoric" multiplicity of the states is auto- 
matically taken into account if the Metropolis criterion 
is tested with the free energy. Furthermore, if all Ni 
are large, the factorials in (||) can be approximated by 
Stirling's formula: 



M 



S{{N}) ^ -fcB^^^nan 
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(8) 



This is the shell-discretized equivalent of the entropy con- 
tribution entering into dJ). By performing this Markov 
process one can thus sample the ion distribution function 
corresponding to the free energy functional (nl) (after an 
initial equilibration time). 

Some final remarks: 

1. The proposed technique is applicable to all geome- 
tries in which the electrostatic energy can be com- 
puted efficiently once the charge density in the vol- 
ume elements is known. In particular, this applies 
to the spherical and cylindrical cell model as well 
as the charged plane, because their high symmetry 
permits equations like (m. 

2. Still referring to particles within a density func- 
tional theory might seem artificial: one could 
equally well transfer small bits of the density be- 
tween shells (under conservation of the integral) 
and employ exp{— /3_F} in the Metropolis criterion 
with the entropy contribution calculated from (g). 
However, Stirling's approximation is only good for 
large Ni, therefore, the global extensive prefactor 
in the exponential function must be large enough 
- for otherwise one even fails to reproduce a con- 
stant density in the limit of zero charge. Since in 
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FIG. 2. Poisson-Boltzmann charge 
fractions P{r) (solid lines) from equa- 
tion (pi) for a valence mixture system 
with R/ro — 100, A = | eo/ro and 
^s/ro = 1- The lowest curve corre- 
sponds to 100% monovalent counteri- 
ons, while for the uppermost curve all 
ions are trivalent. From bottom to top 
the monovalent ions are gradually re- 
placed by trivalent ones in steps of 10% 
of the charge. The dashed line indicates 
the locus of inflection points in P{r) as 
a function of ln(r) for the PB theory 
with only one species of ions (in this 
case for varying ^), whereas the heavy 
dots mark the actual inflection points 
in the distribution functions. 



the spherical geometry the size of this prefactor is 
not at ones disposition (the total number of coun- 
terions per colloid is proportional to its finite and 
possibly small charge), equation (||) is to be pre- 
ferred to equation (0). 

3. Upon measuring the electrostatic interactions via 
the Bjerrum length, the explicit temperature de- 
pendence drops out of the Metropolis criterion. 

4. Further contributions F' to the free energy func- 
tional (e.g., excluded volume or correlation cor- 
rections) can be included by adding corresponding 
terms —(3AF' in the exponential function of (|7|). 

5. The approach bears some resemblance with the 
method of finite elements, for which the transfor- 
mation from a differential equation to |a_functional 
minimization problem is a central ideaO. 



III. EXEMPLARY APPLICATIONS 

In order to demonstrate the utility of the described tech- 
nique, this section presents two possible applications: 
Within the cell model described above the ionic distri- 
bution functions are computed for (i) a system with ions 
having different valences and (ii) a system with addi- 
tional hard core interactions. 



A. Valence mixtures 

The PB equation in the cylindrical cell model cannot be 
solved analytically for a system containing a mixture of 
ions having different valences. For the MC approach, 
however, this is no problem, since it extends in an ob- 
vious way to such a situation: Each species of valence 



Vj is represented by its own density nj{r) and a corre- 
sponding histogram. MC moves only act within the same 
histogram and the electrostatic energy is computed from 
the total charge density p(r) = eo X], "j"'j('')- Figure || 
shows successive stages of an ion exchange process: The 
fraction of charge within a distance r from the rod axis. 



1 r - 

P{r) = — dr27rrp(r). 



(9) 



is plotted for a system with R/ro = 100, A = | eo/ro 
and ^b/^o = 1, in which the neutralizing monovalent 
counterions are gradually (i.e., in steps of 10% of the 
charge) replaced by trivalent ones - until 100% of the 
charge is carried by the latter. The cell was subdivided 
into 500 shells with boundaries equidistant in Inr, 10® 
MC moves were performed to equilibrate the system and 
on average 5 x 10^ further moves were used for sampling 
the distribution. The number of particles which were in 
minority ranged between 200 and 700. 

As expected, the condensed fraction gets larger with 
increasing fraction of trivalent ions, i.e., the distribu- 
tion functions are shifted upwards. Interestingly though, 
there is also a more subtle effect connected with the 
screening of the rod: It has recently been shown that 
within PB theory the point of inflection in P{r) as a func- 
tion of ln(r) localizes the condensation radius and the 
corresponding condensed fractionB. If only one species 
of counterions is present, the locus of inflection points in 
P{r) (parameterized by ^) is solely determined by the cell 
geometry. In Fig. g it can be seen that for the mixtures 
the actual inflection points are shifted towards smaller 
radii, i.e., the rod charge is screened more efflciently. This 
effect must be attributed to global rearrangements of the 
counterions; the distribution of one species clearly de- 
pends on the distribution of the other, different valences 
are correlated - even on the PB level. A more detailed 
analysis in fact reveals that highvalent ions will gather in 
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FIG. 3. Counterion fraction P{r) 
for a cylindrical cell model with 
R/ro = 100, A = 2eo/ro, Is/ro = 1 
and V = 3. The diameter d of 
the ions (implemented via ([L0|) with 
Umax = 3/27rd^) has been varied 
from to 10 ro in steps of ro, which 
shifts the initial rise in P{r) to 
larger values of r. The straight 
dashed line is a fit to the locus of in- 
flection points relevant for the con- 
densation argument. 



the vicinity of the rod at the expense of lowvalent ones. 
This provides a way of further decreasing the free en- 
ergy, therefore, mixtures can screen more efficiently (i.e., 
within shorter distances) than solutions of only one va- 
lence. Notice finally that the observed shift implies that 
there is no "effective" valence V(.g such that a fictitious 
solution of counterions with valence fcs gives the same 
distribution function as the actual valence mixture. 



B. Excluded volume 

The free energy functional ([|) of the plain Poisson- 
Boltzmann equation assumes no interaction beyond elec- 
trostatics. In particular, it neglects all effects resulting 
from an ionic hard core, which range from a suppression 
of high densities up to a packing of particles. A simple 
way for incorporating at least the density limitation is 
the free-volume approximation, characterized by the fol- 
lowing contribution to the free energy density: 



PfFvin) = -n In (l - n/n^^s^^), 



(10) 



where n,„ax is the maximally allowed density. (A related 
approach has recently been used to iiicorporate steric re- 
pulsion in the PB equation, see Ref.Q.) After including 
( [lO| ) into (y) , an MC simulation for a cell model has been 
performed with R/r^ — 100, A = 2eo/ro, ^b/j'o = 1 and 
V = 3. As the limiting density rimax ~ 3/2Trd^ has been 
chosen (with d/ro G {0, 1, ... , 10}), which reproduces the 
correct second virial coefficient for hard spheres of diam- 
eter d. 500 shells and 2000 particles have been used for 
d = To, while for systems with larger ions the number 
of shells has been reduced and the number of particles 
increased (50 shells and 50000 particles for d = 10 rp) 
in order to avoid discretization effects occurring in the 
small innermost shells. After 200-500 MC steps per par- 
ticle for equihbration, a further 5000-10 000 per particle 



have been used to sample the distributions. The effects 
on P{r) are shown in Fig. 0. 

The distribution function for d — ro deviates only 
slightly from the one with d = 0. The reason for this 
is that the contact density in the absence of a hard core, 
"(''o) — 0.447 rQ , is of the same order as the maximum 



density rii. 



0.477 Tq. However, larger values 



of d bring about the expected changes: the contact den- 
sity is reduced in the vicinity of the rod since the ions 
are pushed outwards. Therefore the sharp initial rise of 
P{r) shifts to larger values of r and softens out. Surpris- 
ingly though, the amount-ai condensation (for d = one 
expects 1 — 1/^v w 83%ll3) is mostly unaffected by the 
pronounced changes in the ion distribution at small r. 
This can be seen from the straight line fit to the locus of 
relevant inflection points: Although the condensed layer 
considerably increases in size (from 14.3 ro up to 52.1 tq), 
the thus quantified condensed fraction decreases by less 
than one percent. 



IV. CONCLUSIONS 

A method for finding the equilibrium ion distribution 
belonging to a free energy functional of the Poisson- 
Boltzmann type has been proposed. It avoids any ref- 
erence to the differential equation corresponding to the 
variational problem, which makes it very easy to adapt to 
new free energy functionals. Its utility has been demon- 
strated in two examples, but its range of applicability is 
clearly much larger: It can also be used in the presence of 
salt, with additional free energy contributions (like, e.g., 
correlation correctiongj) and in different geometries. In 
fact, the basic idea could be useful in the general context 
of functional minimization. 
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